Let us consider the dataframe wages, which containes information about 3294 USA working individuals. The data are taken from the National Longitudinal Survey and are related to 1987. The variable as are listed below and the output of the str command is given

A data frame with 3294 observations on the following 4 variables.
exper
   experience in years
male
  1 male, 0 female
school
   years of schooling
wage
   wage (in 1980$) per hour
region
   Center, North, South
## 'data.frame':    3296 obs. of  5 variables:
##  $ exper : int  9 12 11 9 8 9 8 10 12 7 ...
##  $ male  : Factor w/ 2 levels "female","male": 1 1 1 1 1 1 1 1 1 1 ...
##  $ school: int  13 12 11 14 14 14 12 12 10 12 ...
##  $ wage  : num  6.32 5.48 3.64 4.59 2.42 ...
##  $ region: Factor w/ 3 levels "Center","North",..: 1 1 3 1 1 1 1 1 1 3 ...

The aim of the study is to analyze the potential relationship between the response variable wage and the explanatory variables considered in the dataframe. Describe how to perform a preliminary data analysis on this dataframe, using suitable R commands. Moreover, consider the following plots and discuss the possibility of measuring the variable wage in the logarithmic scale

Distribuzione delle variabili

source("./../functions.R", local = knitr::knit_global())
summary(wages)
##      exper            male          school           wage          
##  Min.   : 1.000   female:1569   Min.   : 3.00   Min.   :  0.07656  
##  1st Qu.: 7.000   male  :1727   1st Qu.:11.00   1st Qu.:  3.62673  
##  Median : 8.000                 Median :12.00   Median :  5.20578  
##  Mean   : 8.042                 Mean   :11.63   Mean   :  5.81639  
##  3rd Qu.: 9.000                 3rd Qu.:12.00   3rd Qu.:  7.30723  
##  Max.   :18.000                 Max.   :16.00   Max.   :112.79193  
##     region    
##  Center:1164  
##  North :1105  
##  South :1027  
##               
##               
## 
numeric.graph(wages, c(4,1,3), qqplot = 4, no.density = c(1,3), draw.normal.curve = c(1,3))

## [1] "Skewness indexes"
##       wage      exper     school 
##  9.5491015  0.1555681 -0.3934935 
## [1] "Kurtosis indexes - 3"
##        wage       exper      school 
## 212.9197746   0.6528604   1.2394911
categorical.graph(wages, c(2,5))

La variabile risposta presenta una elevatissima simmetria e curtosi, con i valori praticamente tutti concentrati a sinistra del grafico con la presenza di due outlier elevatissimi. Il qqplot conferma la presenza di una coda destra molto pesante, mentre i quantili centrali seguono abbastanza bene la distribuzione normale. Si consiglia perciò o di rimuovere gli outlier più lontani oppure di introdurre una trasformata.

exper invece sembra seguire piuttosto bene la curva della distribuzione normale (buona simmetria, curtosi leggermente ipernormale).

Considerazioni simili valgono anche per school (anche qui buona simmetria), la distribuzione di presenta con una frequenza del valore centrale (11) molto elevata e una coda sinistra leggermente appesantita.

Vista la natura di queste due variabili non ha molto senso considerare gli outlier.

wages$logwage <- log(wages$wage)
numeric.graph(wages, 6, qqplot = 6)

## [1] "Skewness indexes"
## [1] -1.11844
## [1] "Kurtosis indexes - 3"
## [1] 4.28369

La trasformata logaritmica migliora notevolmente sia simmetria che curtosi della variabile, tuttavia appesantisce molto la coda sinistra (infatti la campana rimane fortemente ipernormale), per cui eventuali assunzioni di normalità sulla variabile vanno fatte con cautela.

Relazioni tra variabili

boxplots.graph(6, c(1:3,5), wages)

Dai grafici non sembrano esserci delle particolari correlazioni tra log(wage) e le altre variabili.

boxplots.graph(1, 2, wages)

mean(wages$exper[wages$male == 'male'])
## [1] 8.323104
mean(wages$exper[wages$male == 'female'])
## [1] 7.732314

Dal grafico e dal confronto delle medie sembra esserci una (debole) correlazione tra exper e male.

Introduzione di trasformate rimuovendo gli outlier

#Rimuovo outlier
newWage <- wages[wages$wage < 30,]
newWage$logwage <- log(newWage$wage)
newWage$wageT <- transformResponseWithBoxCox(newWage$wage~1, newWage, 4)

## [1] "Ideal lambda for wage : 0.4"
numeric.graph(newWage, c(6,7), qqplot = c(6,7))

## [1] "Skewness indexes"
##     logwage       wageT 
## -1.23982911  0.04971715 
## [1] "Kurtosis indexes - 3"
##   logwage     wageT 
## 4.1493448 0.8880249

La trasformata di Box Cox con \(\labda=0.4\) sembra preferibile rispetto a quella logaritmica.

In order to describe the effect of the factor male on the response log(wage) we may analyze this plot, where the probability distribution of log(wage) is represented by considering the kernel density estimates conditioned on the two levels (1 male, 0 female) of the variable male

Il grafico conferma l’assenza di correlazione tra log(wage) e male, infatti le due curve di densità sono praticamente sovrapponibili.

With the commands mod.0<-lm(log(wage) ∼ male,data=wages) and mod.1<-lm(log(wage) ∼ exper*male, data=wages), two regression models are defined for describing the potential effect of male and exper on the response log(wage). Comment the model fitting outcomes given by the function summary (Hint: consider the fact that the average years of experience in the sample is lower for women than for men).

mod.0 <- lm(log(wage) ~ male,data=wages)
summary(mod.0)
## 
## Call:
## lm(formula = log(wage) ~ male, data = wages)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.0445 -0.3073  0.0544  0.3839  3.0325 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.47475    0.01559   94.59   <2e-16 ***
## malemale     0.21826    0.02154   10.13   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6176 on 3294 degrees of freedom
## Multiple R-squared:  0.03023,    Adjusted R-squared:  0.02994 
## F-statistic: 102.7 on 1 and 3294 DF,  p-value: < 2.2e-16
mod.1 <- lm(log(wage) ~ exper*male, data=wages)
summary(mod.1)
## 
## Call:
## lm(formula = log(wage) ~ exper * male, data = wages)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.0906 -0.3050  0.0560  0.3792  3.0468 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     1.193551   0.057470  20.768  < 2e-16 ***
## exper           0.036367   0.007156   5.082 3.94e-07 ***
## malemale        0.463707   0.079062   5.865 4.93e-09 ***
## exper:malemale -0.032071   0.009518  -3.369 0.000762 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6153 on 3292 degrees of freedom
## Multiple R-squared:  0.03792,    Adjusted R-squared:  0.03704 
## F-statistic: 43.25 on 3 and 3292 DF,  p-value: < 2.2e-16

Il confronto tra i due modelli mostra come l’interazione tra exper e male (rappresentato dalla variabile dummy malemale) sia rilevante nel modello, di conseguenza lo sono anche i loro main effects.

Finally, a complete regression model is fitted mod.2 <- lm(log(wage) ∼., data=wages) and the following output is obtained by the R function summary.

summary(lm(log(wage) ~. -logwage, data=wages))
## 
## Call:
## lm(formula = log(wage) ~ . - logwage, data = wages)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.0008 -0.2821  0.0468  0.3673  3.2337 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.279709   0.090321  -3.097  0.00197 ** 
## exper        0.034618   0.004549   7.610 3.55e-14 ***
## malemale     0.246474   0.020607  11.961  < 2e-16 ***
## school       0.122909   0.006278  19.578  < 2e-16 ***
## regionNorth  0.051107   0.024505   2.086  0.03709 *  
## regionSouth  0.047168   0.024969   1.889  0.05898 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5831 on 3290 degrees of freedom
## Multiple R-squared:  0.1364, Adjusted R-squared:  0.1351 
## F-statistic: 103.9 on 5 and 3290 DF,  p-value: < 2.2e-16

Describe how to interpret these results, and then suggest how to proceed with further analyses.

Il modello qui presentato mostra come tutti le variabili siano regressori rilevanti per logwage.

plot(lm(log(wage) ~. -logwage, data = wages))

Il modello presenta alcune criticità dal punto di vista della normalità dei residui. Le ipotesi di omoschedasticità e linearità sono rispettate.

Modello alternativo

Cosa otteniamo utilizzando sempre la trasformata logaritmica senza outliers

fitLog <- lm(logwage~exper+male+school+region, data = newWage)
summary(fitLog)
## 
## Call:
## lm(formula = logwage ~ exper + male + school + region, data = newWage)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.9969 -0.2792  0.0490  0.3710  1.8645 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.299749   0.089186  -3.361 0.000786 ***
## exper        0.036006   0.004493   8.013 1.54e-15 ***
## malemale     0.243299   0.020352  11.954  < 2e-16 ***
## school       0.123402   0.006196  19.915  < 2e-16 ***
## regionNorth  0.056040   0.024193   2.316 0.020597 *  
## regionSouth  0.045451   0.024669   1.842 0.065500 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5754 on 3285 degrees of freedom
## Multiple R-squared:  0.1405, Adjusted R-squared:  0.1392 
## F-statistic: 107.4 on 5 and 3285 DF,  p-value: < 2.2e-16
plot(fitLog)

drop1(lm(wageT~exper+male+school+region, data = newWage), test = "F")
## Single term deletions
## 
## Model:
## wageT ~ exper + male + school + region
##        Df Sum of Sq    RSS    AIC  F value    Pr(>F)    
## <none>              3300.7  21.66                       
## exper   1     55.62 3356.3  74.66  55.3587 1.275e-13 ***
## male    1    170.70 3471.4 185.60 169.8906 < 2.2e-16 ***
## school  1    444.86 3745.5 435.76 442.7430 < 2.2e-16 ***
## region  2      3.29 3304.0  20.94   1.6363    0.1949    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Dal risultato di drop1 risulta che per la nuova trasformata il regressore region non è significativo.

fitBoxCox <- lm(wageT~exper+male+school, data=newWage)
summary(fitBoxCox)
## 
## Call:
## lm(formula = wageT ~ exper + male + school, data = newWage)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.7995 -0.6243 -0.0142  0.6259  4.4220 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.997868   0.153070  -6.519 8.16e-11 ***
## exper        0.057779   0.007824   7.384 1.93e-13 ***
## malemale     0.462388   0.035457  13.041  < 2e-16 ***
## school       0.226935   0.010796  21.021  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.003 on 3287 degrees of freedom
## Multiple R-squared:  0.1523, Adjusted R-squared:  0.1516 
## F-statistic: 196.9 on 3 and 3287 DF,  p-value: < 2.2e-16
plot(fitBoxCox)

La trasformata calcolata con il metodo di Box-Cox permette di rientrare (seppur parzialmente) nell’ipotesi di normalità dei residui.

AIC(fitLog, fitBoxCox)
##           df      AIC
## fitLog     7 5710.026
## fitBoxCox  5 9362.390
BIC(fitLog, fitBoxCox)
##           df      BIC
## fitLog     7 5752.719
## fitBoxCox  5 9392.884

Gli indici \(R^2\) sono più alti nella trasformata Box Cox rispetto al logaritmo, tuttavia le statistiche AIC e BIC tendono a preferire quest’ultima.